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‘Aerospace Engineer, NASA Goddard Space Flight Center, Greenbelt, MD, 20771. 


1 


NOMENCLATURE 


-^/C5 


r, u, x 
w(t) 


m 

p(ti,tj) 

Q(t),Qd{t, u) 

P(r,t) 

V* 


Pc{t),Pnc{t) 


Pc(^> U ) , Pnc{ ^i) 


S(t) 


Position, velocity, and position/velocity state vector of spacecraft 
k = 1,2 relative to central body 

Position, velocity, and position/ velocity state vector of spacecraft 2 
relative to spacecraft 1, r = R 2 — jRi, etc. 

Model of the unperturbed time evolution of x, x(t) = /(x,t) 
Hypothetical noise-like process perturbing the time evolution of x 
according to x{t) = /(x, t) + B{t)w{t) 

Error in the relative state, e = x — E[x] 

Relative state error covariance matrix, P(t ) = E[e(t)e(t) T ] 

Relative state error auto-covariance matrix, P(U, tj ) = E [e(ti)e(tj) T ] 
State transition matrix for e, 4>(t, t{) = A(t)$(t,ti), A(t) = 

d f /dx,${ti,ti) = I 

Power spectral density of wit), and process noise covariance, 
Qd{t,U) — ft { $(u,ti)B(u)E[w(u)w(s) T ]B(s) T $(s,ti) T dsdu 
Probability density function of the relative position 
Convex set of points in R 3 representing a protection volume sur- 
rounding the combined hard-body volumes of the two spacecraft, at 
a particular time, U 

Instantaneous probability of collision, p c (U) = Pr (r(ti) G V*), and 
probability of not colliding, p n c(U ) — Pr(r*(£i) 0 Vi) 

Cumulative probability of collision, and probability of not colliding, 
over a time interval 

Auto-covariance of joint density associated with p n c(t,U) 

Survivor function, s(t) — Pr(£ c > £), where t c is the time of collision; 
is equal to p nc (t, U) for ti = 0. 


INTRODUCTION 

The problem of estimating the probability of collision between an operational spacecraft 
and orbital debris has been widely studied. Foster and Estes 1 and Akella, et al. 2 describe 
the debris problem, and Chan’s approximate analytical solution, 3 which is valid for cases 
in which the uncertainty in the relative positions is much larger than the combined hard- 
body radius, is now widely accepted. All of these methods assume that the encounter is 
brief, so that the relative motion may be assumed to be along a straight line normal to the 
encounter plane, and focus on the computation of the instantaneous probability of collision 
at the single moment of closest approach. Patera 4 describes one way to relax the “linear” 
motion assumption, and also accommodate time- varying uncertainty. 

An implicit, but not necessarily obvious, assumption in all of these techniques is that the 
relative position uncertainty is perfectly correlated in time. If there is any mis-modeling of 
the dynamics in the propagation of the relative position error covariance matrix, time- wise 
de-correlation of the uncertainty will increase the probability of collision over a given time 
interval. Some examples given in the sequel illustrate this point. If the dynamics uncertainty 
can be accurately modeled as a Brownian motion diffusion process, the time-evolution of the 
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probability density of the relative position can be modeled by the Fokker-Planck Equation 
(FPE), 

The FPE however is very difficult to solve efficiently for systems of high dimensionality. The 
recent work of Kumar, et al . 5 has proposed the use of more efficient solution techniques, 
but this work has not yet achieved the goal of applying these techniques to a realistic space 
encounter. 


In current practice, therefore, the Monte Carlo method is still widely used to study 
problems associated with low- velocity encounters. This approach is not without issues of 
its own. A question that is open to discussion is always how many Monte Carlo trials 
need to be run to assure adequate results. One way to address this question is to specify a 
confidence interval, based on the number of trials, but once again, the choice of confidence 
limits is left open to the analyst. A more subtle issue is right censoring, which is associated 
with the concept that one never knows if a particular trial would have ended in collision 
if the simulation were run a bit longer. A further issue with the Monte Carlo technique 
is that it is only as good as the accuracy of its inputs; uncertainties in the assumed error 
models will pollute the results. In this technique, the initial conditions are perturbed based 
on an initial covariance matrix. In an operational system, this covariance is typically the 
result of an orbit determination (OD) process. Whether or not the OD itself utilized pro- 
cess noise, the accuracy in predicting the resulting covariance matrix to a future epoch will 
be limited, and empirical techniques for the addition of process noise to these covariance 
propagations have been found to enhance the prediction accuracy . 6 The physical justifi- 
cation for adding process noise is that certain parameters in the dynamics model used for 
state and covariance propagation are not known to high precision. Therefore, the addition 
of diffusive uncertainty, which the process noise accomplishes, introduces a conservative 
bound on the resulting propagation errors, so long as the process noise covariance itself 
is accurate. Setting aside for the moment the issues of finding accurate initial covariance 
and process noise covariance, and whether or not white noise is the best way to model the 
dynamics uncertainty, there remains the issue that Monte Carlo analysis often does not 
introduce additional noise at each simulation step based on the process noise covariance. 
This omission would appear to be questionable, since without it, the propagations in the 
Monte Carlo trials will omit any uncertainty in the dynamics that is actually present. 

This work takes the view that, for the present, Monte Carlo analysis is the best available 
tool for handling low-velocity encounters, and suggests some techniques for addressing the 
issues just described. One proposal is for the use of a non-parametric technique, known 
as the Kaplan-Meier estimator, that is widely used in the realm of actuarial and medical 
studies for handling the right censoring bias inherent in Monte Carlo trials. The other 
suggestion is that accurate process noise models be used in the Monte Carlo trials to which 
the non-parametric estimate is applied. A further contribution of this paper is a description 
of how the time-wise de-correlation of uncertainty increases the probability of collision. 

The remainder of this paper is organized as follows. First, a fairly extensive background 
section describes the issues highlighted above, and illustrates the issues with two examples 
of increasing complexity. The next section describes a proposed procedure for forming non- 
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parametric estimates of the collision probability, and applies the procedure to the examples 
of the background section, as well as a somewhat realistic formation flying problem. The 
final section concludes, summarizes, and offers suggestions for further related work. 

BACKGROUND 

This section first derives the joint density associated with Pncit^U) over a series of 
discrete sample points in order to show the effect of time- wise de-correlation on p c (t, £*), 
due to the diffusive nature of process noise. It then illustrates this effect with two simplified 
examples of reduced dimensionality so that the geometry involved is evident. Next, a 
subsection describes some existing methods for determining accurate values for the initial 
error covariance and process noise covariance. The final subsection discusses some concepts 
from the actuarial and life sciences literature known as the hazard and survivor functions, 
and discusses their relevance to the problem of determining the probability of collision over 
an interval. 

Joint Density 

Considering a discrete series of time samples, the probability of collision over an interval 
of such samples, p c {t n , to)? is the probability associated with the union of the events that the 
perturbed relative positions, n, are inside the volumes, V*, defining the protection region 
at each sample time t*. This may be written as the complement of the probability of no 
collision, Pnc(t n ,t o), associated with the intersection of the events that the r\ are outside 
the protection regions, 


Pc(tni to) — Pr(LbeVi) (2) 

i=0 

n 

= l - Pr(p| n e Vi) (3) 

i = 0 

— 1 ~ Pnc(tn)t 0 )) (4) 

where V* represents the complement of the set V*. Note that Eq. 2 is only an approximation, 
and unless a fairly dense sampling of time points are considered, Eq. 2 might underestimate 
Pc(tni to) severely, since it could miss collisions that occur between the sample times. 


With this caveat in mind, the joint probability that approximates the probability of 
not colliding over an interval is given by the integration of the joint density function 
p(r n , r n- l, • • * 7T, ro, E) over each of the complements of the protection regions, 


Pncitn-) to) = dr n dr n -i • • • dr i I dr 0 p(r n ,r n -i, 

J J Jy n J J JVn-i J J Jv 1 J J Jy o 


•ri,r 0 ,E), 

( 5 ) 


where, assuming a Gaussian density, p(r n , r n _ i, • • • ri, ro, S) is 


pi^ni r n— 1? ' * ‘ ?T? ^0? S) 


exp (-^ [el, e^ 1; • • ■ ef, eg'] E 1 [e n , e n _ x , e 0 ]) 

V / (27r) 3 ("+ 1 )|E| 


(6) 
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and where ei = ri — E [r*] are the errors in the relative positions. The autocovariance of the 
relative position errors, E, is 


E = 


Pr n 

Pr n r n - 1 

Pr n ri 

Pr n ro 

Prn-ir n 

Prn-, 

Pr n -iri 

Pr n -ir 0 

Prrrn 

p 

A r\r n -i 

Pr i 

P riro 

Pr 0 r n 

Pror n - 1 ‘ ' 

P ro n 

Pro 


(7) 


where the covariance of the relative position error generally includes a cross-covariance 
between the the errors in the satellites’ absolute positions, 


Pri — PR\ (j'i) H"-Pr 2 (^) PR\R2 i^i) Pr\R2 (^) 

and the relative position autocovariance over a time interval is 


(8) 


tj ) Pj 

i 

$r(ti,to)Po$r(tj,to) + E* r tk)Qd(j'ki ^k— l)* 

k=j 


(9) 


In Eq. 9, denotes the upper three rows of the 6x6 state transition matrix, 

i.e., Using Eq. 9, one may decompose the expression for E into a 

summation: 



to) 

&r(tn—li t 0 ) 


*&r(tni to) 
(tn— 1) *o) 

T 



*&r(tni t±) 

E = 

$r(t l,io) 
I 

Po 

$r(t l, i 0 ) 

7 

+ 

$r(t2,h) 

I 

0 

Qd(tl,to) 

&r(t2,t\) 

I 

0 



*&r(tm tn— l) 


^r( 7 n» tn— l) 

T 

' 7 ' 

r 


7 


7 


0 


+ ••• + 

0 

Qd(tn— lj ^n— 2 ) 

0 

+ 

0 

Qdi^n ) ^n—l) 


0 


0 


_ 0 _ 



-i T 


(10) 


Clearly, the evaluation of Eq. 5 poses challenges perhaps as great as the solution of the 
FPE, and at best, will only give an approximate result. It does however offer insight into 
the nature of the problem. Since the first term in Equation 10 has at most the rank of Po ? 
and each of the remaining terms have at most the rank of their corresponding Qd(U,U- 1 ), 
then it is clear that |E| = 0 if there is no process noise. This merely indicates that, given 
a realization of the random initial relative position, r*o, the r* are not random variables 
when there is no process noise. In this case, the initial distribution of the position errors 
completely determines p n c(tn,t o) — Pr(p|£_ 0 n ^ V$). Thus, when there is no process 
noise, the joint probability that there are no collisions over a series of discrete sample times 
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reduces to the individual probability that the perturbed initial relative position does not lie 
within the region representing the combined mapping back to the initial time of the various 
protection regions at each sample time, that is, 

n n 

Prid^m to) = Pr(f| n e Vi) = Pr(p| r 0 € g(V 0 ,V i, ...,%)) (11) 

i = 0 i=0 

where g is the mapping of the protection regions to the initial time. Note that this mapping 
should be performed with the dynamics governing the PDF, which are not necessarily the 
same as the dynamics governing the time evolution of the protection regions (although 
typically they are the same for the satellite collision problem). 

Effect of Time Correlation 

Equation 10 illustrates two sources for time-wise de-correlation to arise. The first term 
in the summation shows that the uncertainty in initial conditions will evolve according to 
the system dynamics, which will induce a time-wise de-correlation on the order the time 
constant of the dynamics. The remaining terms in the summation show how the presence 
of process noise induces a persistent de-correlating effect. Figures 1 and 2 illustrate these 
points with a pair of examples. To simplify the illustration, these examples do not use 
realistic two-body orbital dynamics; a later example introduces more realistic two-body 
motion. From these examples, it is clear that time correlation can drastically affect the 
outcome, even when the instantaneous moments of the density are identical. 

In the example Figure 1 illustrates, two spacecraft in circular orbits are initially sepa- 
rated by 2.5 deg in true anomaly, and drifting together at a rate of ldeg per orbit period. 
This example only considers the along-track motion of the two spacecraft. At each point 
in time, there is a unit variance Gaussian distribution of the relative position errors. In 
each of the three subplots, a different time correlation holds. A ±0.5 deg region represents 
the projection onto the approach path of a substantial protection region in excess of the 
combined hard-body radius, which fills a parallelogram when time is the horizontal axis, 
as in this figure. The dark circles impinging the protection region represent close approach 
violations. 

In the example Figure 2 illustrates, a spacecraft orbits a neighboring circular orbit 
spacecraft in a plane inclined 30 deg to the local horizontal so that its relative motion is a 
circle, with radius 2.5 km. This example only considers the two-dimensional motion of the 
spacecraft in the orbital plane. As in the one-dimensional example, at each point in time, 
there is a unit variance Gaussian distribution of the relative position errors, and as before, 
in each of the three subplots, a different time correlation holds. A 1.5 km radius circular 
region represents the “protection region” in excess of the combined hard-body radius, which 
creates a spiral tunnel when time is the third axis, as in this figure. As before, the dark 
circles impinging the protection region represent close approach violations. 

To further illustrate the effect on collision probability of changing the time correla- 
tion, Table 1 shows some results of directly evaluating the integrals of Eq. 5, for the one- 
dimensional example that Figure 1 shows. In order to make this calculation practical, the 
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Timewise Uncorrelated - Bernoulli Trials 



Figure 1: One-dimensional along-track approach trajectory example (1000 trials, 100 
time increments). 





Figure 2: Two-dimensional circular relative motion example (1000 trials, 100 time in- 
crements) . 


7 


number of time increments had to be reduced from 100 to 10 steps 2 * * . The table clearly shows 
that decreasing the time correlation substantially increases the probability of collision over 
the interval, for this example. 


Table Is Effect of changing correlation time (vs. baseline value of 0.3 sec), for the in- 
termediate case from Figure 1, with ten (vs. 100) time increments. 


Correlation Time [sec] 

PnciX') 0) 

Pc(l,0) 

0.01 

0.412 

0.588 

0.03 

0.416 

0.584 

0.1 

0.489 

0.521 

0.3 

0.611 

0.389 

1.0 

0.740 

0.260 

3.0 

0.805 

0.195 

10 

0.834 

0.166 


Covariance Prediction 

In order for any type of collision prediction to be accurate, there is a need to predict 
the covariance accurately, which requires an accurate process noise covariance. This process 
noise covariance need not be based on a white noise model, e.g. Wright’s work on sequentially 
correlated process noise based on gravity mismodeling; 7 what matters is getting a good 
prediction. In some cases, empirical techniques, e.g. Duncan and Long, 6 are currently in 
use. An apparently little- known fact is that sequential smoothing techniques can provide 
optimal estimates of the process noise sequence. Appendix X of Bierman’s text 8 shows 
the estimate for w(t). Recently, Psiaki described a sigma-point smoother 9 that also shows 
the update for the process noise and its covariance. In the usual sequential form of the 
Rausch-Tung-Striebel smoother, the process noise update takes the following form: 


j. _ p.fo T p- 1 

(12) 

Jiw = QdiPi+xi 

(13) 


(14) 


(15) 

p* = Pi- Jix{p* + i - Pi+i)JL 

(16) 

Qdi = Qdi ~ Jiw(P t +\ Pi+l) Jiun 

(17) 


where (•), and (•) denote the specified quantities just prior to, and after, the forward filter 
measurement update, respectively, and (•)* denotes the backward smoother estimate. One 
should note that this technique does not directly provide Qdi, the forward filter process noise 
covariance, which is the parameter needed for accurate covariance prediction. However, the 
estimates it provides of the process noise sequence, u;*, may prove useful as part of the filter 
tuning process. 

2 One consequence of this large time step is that, for the smaller correlation times, it becomes increasingly 

likely for the Gauss-Mar kov model to take steps completely across the protection region between time steps; 

such steps are considered to be collisions. 
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Actuarial Analysis of Failure Times 

A common problem in the fields of actuarial and medical studies is the determination 
of models for the time until some event, most often death, occurs. In these studies, there 
is typically a set of observations of the time of this event for a large group of individuals, 
and the analyst wishes to determine a distribution for these data. These data typically 
do not cover the entire lifespan of the individuals, so that many of the individuals are still 
surviving at the end of the data span. It may also occur that some individuals leave the 
study early, for reasons other than death. The literature refers to these conditions as right 
censoring. There are related concepts of left censoring and left truncation. In left censoring 
an individual failure occurs prior to some specific time, but the exact time of failure is 
unknown. In left truncation, an individual that would otherwise have been included in the 
study is not because they have already died prior to the beginning of the data span. 

According to a standard textbook, 10 there are three specifications for the failure time 
distribution that are “particularly useful:” the survivor function, the probability density, 
and the hazard function. The survivor function, s(t), is the probability that the time of 
the event in question is greater than some specified value, so it is the complement of the 
cumulative distribution function. The survivor function is particularly important because 
“the right tail [of the distribution] is the important component for the incorporation of 
right censoring.” The hazard function is the “instantaneous rate at which failures occur for 
items that are surviving” at a given time, for continuous time models, and the “conditional 
probability of failure at [a specified time], given that the individual has survived to [the 
specified time]” for discrete-time models. General models for survival time can have both 
continuous and discrete components to the hazard function. 

The optimal non-parametric estimate for the survivor function is found by maximizing a 
likelihood function on the space of all survivor functions. This is a discrete survivor function 
formed as the product of hazard components, 1 — hi/r^ where hi is the number of items 
that fail at each and r{ is the number of items remaining at risk for failure at £*. Note 
that this assumes that the lifecycles of individuals are independent from each other, and 
that any censoring that occurs is independent, but does not assume independence between 
between time samples for a given individual, which the Bernoulli trials assumption often 
used in engineering failure analysis imposes. The result is known as the Kaplan-Meier or 
product limit estimator for the survivor function: 




( 18 ) 


Note that if there is no right censoring at any interior time samples (i.e. no one leaves the 
study prior to its completion), then r*+i — r{ — hi so that the survivor function reduces to 
the intuitive value of s(tj ) = 1 — Yjl=i hi/m, where m is the total number of individuals in 
the data set. 


Applying the obvious relationships of these concepts to the satellite collision problem is 
the subject of the next section. 
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NON-PARAMETRIC ESTIMATE FOR p c (t n ,t 0 ) 

This section proposes a procedure based on the discussion above for computing collision 
probability over an interval, that addresses all of the issues mentioned previously. Next, it 
applies the procedure to the examples of the previous section, and somewhat more realistic 
example of a formation flying problem. 

Proposed Procedure 

Supposing that the Kaplan-Meier survivor estimate resulting from an accurately per- 
formed Monte Carlo study is a suitable proxy for the probability of not colliding over the 
interval of the Monte Carlo study, 

Pc(t 715 *o) = 1 Vndj'm *o) = 1 Pr(£ c > t n ) = 1 J(t n ), (19) 

this paper suggests the following procedure for estimating the probability of not colliding 
over that interval, and hence an estimate of p c (£ n ,£ o). 


1. Obtain realistic estimates of both the initial covariance and the process noise covari- 
ance associated with the prediction model, via empirical techniques such as Duncan 
and Long 6 suggest, or via a smoother, such as the sigma-point smoother described by 
Psiaki. 9 


2. Perform a monte-carlo prediction study with m trials, using these realistic values for 
both the initial covariance and the process noise. 


3. Compute the Kaplan-Meier estimate of the survivor function, 11 i.e. the probability 
that the time of collision is later than a given time: 


%) = n 

2=1 


Hi 


(20) 


where r* is number of trials at risk, i.e. those trials that have not resulted in a collision 
up to time tj, (i.e. those trials that would be right censored if the trial ended at £j), 
and hi is the number of trials that have resulted in a collision at time tj. Note that for 
most Monte Carlo studies, there should be no interior right censoring, so the survivor 
function will reduce to 

3 

s(tj ) = 1 — ^ hi/m. (21) 

2=1 


4. Use Greenwood’s formula 12 to compute the sample variance of s : 


var[5(^)] = s(t J ) 2 ^ — 


“ ri{ri - hi) ’ 


(22) 


5. If s is near one (or zero), so that confidence limits based on the Greenwood sample 
variance exceed one (or zero) , then use instead the technique suggested by Kalbfleish 
and Prentice: 10 compute 

i 2 (t) = vai'[log s(»]/ log[s(i)] 2 , (23) 
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which is the asymptotic variance of log[— log s(£)], and then for the confidence interval 
of s(t) use 

[s(t)] exp(± *^ W , (24) 

where z a is the normal deviate corresponding to the desired confidence interval, e.g. 
z a = 1.96 for 1 — a = 95% confidence limits. 

6. If s is exactly equal to one, then more trials are needed, since by conventional wisdom 
no system has immortal life. If running more trials is not practical, then an estimate, 
which is based on the binomial distribution, of the reliability that would correspond 
to 1 — a two-sided confidence limits, with the right confidence interval having a/2 of 
the mass resting solely at 100% reliability, is given by 

Pnc = exp(log(a/2 )/m), (25) 

(the left confidence interval may be found from the binomial cumulative density eval- 
uated at a/2, or via the normal approximation to the binomial density via the usual 
method.) 

7. If the confidence limits are larger than desired, run more trials. Note that the equation 
above may be solved for m to specify the number of trials required to establish a given 
reliability, with given confidence bounds, if the given reliability is close to one. 


Results for Constant Density Examples 

Figures 5 and 6 illustrate the results of applying the procedure outlined above for the 
simplified examples that Figures 1 and 2 show. In addition, Figure 7 shows the results 
of applying the procedure to a simple three-dimensional example (circular relative motion 
trajectory test case, with R = a = 1 and r = 0.3) that Reference 4 describes. 

For the one- and two-dimensional perfect time correlation cases, it is possible to compute 
the exact value of p c {t n , to), since the mapping of the protection region described in Eq. 11 is 
simply the projection of these regions onto the t = 0 plane; in other words, it is the shadow 
cast by the protection regions back onto t = 0, as Figures 3 and 4 show. In Figure 3, the 
shadow of the two-dimensional protection region onto t — 0 is one-dimensional, and the 
third dimension shows the constant probability density. In Figure 4, the shadow of the 
three-dimensional protection region onto t = 0 is a two-dimensional annulus, and it is not 
possible to show the probability density, as this would require a fourth dimension. These 
figures illustrate the point made by Eq. 11 and its accompanying text. 

For these cases, the shadows cast especially simple geometries. In Figure 3 the shadow 
projects onto one of the tails of the normal density that governs the initial error distribution, 
so that computing p c reduces to evaluating the cumulative normal density over the portion 
of the tail outlined in Figure 3. In Figure 4, the shadow is a disk with a missing center, 
so that computing p c reduces to computing the circular error probable, or the cumulative 
Rayleigh distribution over the outer disk, and then subtracting the Rayleigh mass contained 
in the inner disk, or hole. It is convenient to find these probabilities by evaluation of the 
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Figure 3: Shadow cast on t = 0 
plane by motion of protec- 
tion region in Figure 1 



Figure 4: Shadow cast on t = 0 
plane by motion of protec- 
tion region in Figure 2 


X 2 density, with v degrees of freedom, 

f z z^~ 2 ^ 2 e ~ z / 2 

= / jo — ; — d z, (26) 

x V ’ ; J 0 2*V 2 r(i//2) 5 v ' 

since this density generalizes the symmetric Gaussian density to arbitrary dimensionality. 
For the one-dimensional example (is = 1), 

Pnc( 1, 0) = 1 - p c ( 1, 0) = 1 - * y 2 = 0.843, (27) 

where the integral must be halved since only one tail is needed, and for the two-dimensional 
example (is = 2), 


Pnc( 1, 0) = 1 - p c ( 1, 0) = 1 - m x 2 (4 2 , 2) - m x2 (l, 2) = 0.394. (28) 


For the three-dimensional example, Reference 4 gives the following result for the exact 
collision probability for the perfect time correlation case: 


Pc(t,ti ) = - W-exp 
<T V 7T 


p c (i,o) = 


(i) 


exp 


(i? 2 + r 2 )l f r . 

/si 
. Jo 

1/ 

J Jo 


2a 2 
— (I 2 + 0.3 2 )" 1 


sinh 


R\Jr 2 — x 2 




dx, 


2(1)2 


sinh 


/ (l)V0.3 2 -x 2 

V 


dx = 0.066146, 


Pnc(l, 0) = 1 — p c (l, 0) = 0.934. 


(29) 

(30) 

(31) 


The figures show that the proposed procedure estimates the known probability for the 
perfect correlation cases to within the confidence limits. For intermediate and zero correla- 
tions, the procedure produces the expected reductions in the survival probability. 
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Survival Probability - Zero Time Correlation 



Figure 5: The empirical survival probability (blue), with ±95% confidence interval (red 
and green), for the one-dimensional case Figure 1 shows. For the perfect time 
correlation case, the exact value after one orbit period is 0.843. 


Survival Probability - Zero Time Correlation 



Figure 6: The empirical survival probability (blue), with ±95% confidence interval (red 
and green), for the two-dimensional case Figure 2 shows. For the perfect time 
correlation case, the exact value after one orbit period is 0.394. 
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Survival Probability - Zero Time Correlation 



Survival Probability - Intermediate Time Correlation 



Fraction ot Orbit Period 


Figure 7: The empirical survival probability(blue), with ±95% confidence interval (red 
and green), for the three-dimensional case Reference 4 describes (1000 trials, 
100 time increments). For the perfect time correlation case, the exact value 
after one orbit period is 0.9339. 
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Results for a Realistic Formation Flying Example 

Figures 8 and 9 shows a somewhat realistic three-dimensional formation flying example. 
A spacecraft orbits a neighboring circular orbit spacecraft in the neighbor’s orbit plane, so 
that its relative motion is a 2 x 1 ellipse, with semi-major axis 200 m. There is uncertainty 
in both the initial relative position and velocity (1 m and 1 mm/sec, lcr, respectively, with 
— .95 correlation coefficient between radius and speed, and between radial velocity and in- 
track position), as well as acceleration process noise of 3 milligee/sec 1 / 2 . The protection 
region is a 50 m sphere. The example is only partially realistic, since it assumes two-body 
gravity, and the process noise is artificial white noise. In reality, the process noise arises due 
to mis-modeling of the dynamics. Figure 9 shows an expected step decrease in the survival 
probability each time the relative trajectory crosses the orbit of the reference spacecraft. 

CONCLUSIONS 

This paper has proposed a procedure for computing the collision probability over an 
interval using a non-parametric technique, known as the Kaplan-Meier estimator, that is 
widely used in the realm of actuarial and medical studies. The estimator accommodates 
the right censoring bias inherent in Monte Carlo trials. This proposal will only be as good 
as the fidelity of the Monte Carlo procedure, and this paper has offered some suggestions 
to improve the accuracy of the process noise and initial covariance matrices used in the 
Monte Carlo trials to which the non-parametric estimate is applied. This paper has also 
described how the time-wise de-correlation of uncertainty due to process noise increases the 
probability of collision. 

Although other techniques have been proposed for low- velocity encounters that are less 
costly to perform than Monte Carlo analysis, unless there is sufficient justification for ignor- 
ing uncertainties in the propagation model, the use of any technique that does not properly 
model time-wise de-correlation due to the model errors would not appear to be advisable. 
One possible avenue for future work that could address this problem might be to generalize 
works such as Reference 4 so as to perform a four-dimensional integration over space and 
time. This would be analogous to integration over the volume of the spiral tube that Fig- 
ure 2 illustrates, but generalized to a hyper-tube through three spatial axes perpendicular 
to the time axis. Alternatively, continued progress in the numerical solution of the FPE 
may offer the best hope for efficiently solving this difficult problem. 

ACKNOWLEDGMENTS 

Discussions with the following individuals have shaped the development of this work: 
Landis Markley, Terry Alfriend, Conrad Schiff, Marco Concha, Mrinal Kumar, John Junk- 
ins, Mark Psiaki, and Carol Moore. 

REFERENCES 

1. J. L. Foster, Jr. and H. S. Estes, “A Parametric Analysis of Orbital Debris Collision 
Probability and Maneuver Rate for Space Vehicles,” Tech. Rep. JSC-25898, NASA 
Johnson Space Center, Houston, TX, 1992. 


15 



Figure 8: Three-dimensional in-plane relative motion example (100 trials). 
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Figure 9: The empirical survival probability (blue), with ±95% confidence interval (red 
and green), for the three-dimensional case Figure 8 shows (100 trials). Each 
of the step decreases in s(t) occurs when the relative trajectory crosses the 
orbit path of the reference spacecraft. 
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